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ABSTRACT 

This  report  describes  the  procedures  used  to  create  0  and  C-type  solution 
grids  around  arbitrary  two-dimensional  bodies.  The  grid  generation  procedure 
is  a  modification  of  the  0-type  hyperbolic  method  of  Steger  and  Chaussee^. 

The  procedure  described  in  Reference  (1)  was  modified  by  including  additional 
dissipation  terms  and  by  changing  the  form  of  the  dissipation  described.  The 
modifications  necessary  to  produce  a  C-type  grid  are  discussed  and  examples 
are  provided. 

A  brief  description  of  the  theory  and  development  of  the  governing 
hyperbolic  equations  is  provided.  Example  results  for  both  an  airfoil  section 
and  a  complex  body  (X-24C)  are  shown.  A  discussion  of  the  user  definable 
input  variables  and  their  effects  on  the  resulting  grid  is  included. 

For  many  solution  procedures  it  is  desirable  that  the  distribution  of 

O 

points  around  the  airfoil  be  "second  order  smooth."  Since  most  airfoils  are 
defined  with  a  somewhat  random  distribution  of  points,  a  routine  was  used  to 
redefine  the  coordinates  in  a  smooth  distribution.  The  routine  also  allows 
selective  clustering  of  grid  points  in  regions  of  interest.  A  brief  description 
of  this  routine  and  the  effects  of  the  different  input  parameters  are  included 
as  Appendix  A. 
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I.  INTRODUCTION 


Finite  difference  solutions  of  the  partial  differential  equations  that 
govern  fluid  flow  about  arbitrary  bodies  require  the  development  of  a  computa¬ 
tional  grid  around  the  body.  In  particular,  a  grid  system  that  conforms  to 
the  shape  of  the  body  has  the  very  important  advantages  of  simplifying  the 
partial  differential  equation  solution  technique  and  simplifying  the  satisfaction 
of  the  boundary  conditions  at  the  body.  In  many  cases,  the  character  of  the 
entire  flow  field  is  determined  in  the  high  gradient  regions  near  the  body. 
Therefore,  the  development  of  a  body-conforming  coordinate  system  is  a  necessary 
and  very  important  first  step  in  numerical  solutions  of  fluid  flow  problems. 

Algebraic  transformation  techniques  for  grid  generation  have  been  used 
for  some  applications  and  have  the  advantage  of  being  relatively  simple  and 
fast.  However,  the  most  popular  techniques  use  the  solution  of  partial 
differential  equations  (PDE).  Procedures  exist  that  solve  elliptic,  hyperbolic 
and  parabolic  PDFs.  Elliptic  routinesare  most  often  used  and  probably  are 
more  generally  applicable  to  a  wider  range  of  problems.  Hyperbolic  routines 
are  relatively  fast  and  simple,  but  are  not  easily  adapted  to  physical  far 
field  constraints  (such  as  internal  flows  and  multiple  bodies).  Parabolic 
routines  have  been  recently  developed  and  show  some  promise.  The  PDE  methods 

attempt  to  produce  nearly  orthogonal  grids  with  assurance  of  no  crossover  of 
adjacent  grid  lines. 

This  report  will  discuss  the  characteristics  of  a  hyperbolic  grid  generation 
routine  modified  to  produce  a  computational  grid  over  arbitrary  two-dimensional 
bodies.  The  hyperbolic  grid  generation  routine  was  chosen  because  it  is  quite 
fast,  provides  nearly  orthogonal  grids  and  has  good  user  control  of  the  grid 
spacings.  The  routine  should  work  well  for  generating  a  grid  around  bodies 
where  there  are  no  physical  constraints  (such  as  other  bodies,  walls,  etc.) 
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within  the  region  of  the  grid.  The  actual  grid  generation  starts  at  the  body 
at  user  defined  locations  and  marches  out  to  an  outer  boundary. 

The  core  of  this  program  was  adapted  from  the  procedure  described  by 
Steger  and  Chaussee  in  Reference  1.  The  purpose  of  this  report  is  to  describe 
and  document  the  modifications  made  to  the  0-grid  algorithm  of  Steger  and 
Chaussee  and  to  provide  examples  of  how  the  various  user  defined  variables 
affect  the  results. 

Another  program  called  SMOOTH  has  been  developed  to  allow  a  smoothly 
varying  distribution  of  the  points  describing  the  body.  The  same  process  that 
insures  smoothness  also  allows  arbitrary  clustering  and/or  spreading-out  of 
points  around  the  body.  This  program  is  briefly  described  and  some  representa¬ 
tive  results  are  provided  in  Appendix  A. 
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II.  THEORY  DEVELOPMENT 


In  the  grid  generation  algorithm  developed  by  Steger  and  Chaussee, 
definitions  of  grid  orthogonality  and  transformation  Jacobian  were  chosen  to 


devise  a  scheme  mapping  (x,y)  to  (^.n). 

A  typical  transformation 

for  an 

airfoil  is  shown  in  Figure  1. 

*5^  *  V.  -  » 

Orthogonal ity 

(1) 

Inverse  Jacobian 

(2) 

Note  that  in  general,  area  integrands  in  physical  and  computational  domains 
can  be  written 


dA  =  dXdY  =  1/J  dcdn 

Numerically  a?  =  An  =1,  so  the  inverse  Jacobian  approximates  the  physical  cell 
area.  Equations  (1)  and  (2)  are  locally  linearized  using 

X  =  X°  +  AX 
Y  =  Y®  +  ay 

where  X°,  y°  is  a  nearby  location.  The  linearized  set  of  equations  become 
(after  some  algebraic  manipulation) 
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COORDINATE  TRANSFORMATION 


FIGURE  I 
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or 


A  +  B  =  f  (3) 

The  specific  numerical  procedure  for  solving  equation  (3)  depends  on  the  class 

of  equations  being  solved  (elliptic,  hyperbolic,  or  parabolic).  The  class  can 
be  determined  by  inspecting  A. 


where 

j  =  X  +  Y 

-1  ^ 

Since  B  A  is  symmetric,  it  has  real  eigenvalues,  specifically: 


^1.2  =  -  +  (x^o  y^o  ^  Y^o 

This  indicates  that  the  system  is  hyperbolic  and  can  be  marched  in  n. 

The  finite  difference  solution  scheme  used  is  centrally  differenced  in  ^ 
and  backward  differenced  in  n.  The  scheme  can  be  written  as 


*  (’j 


Where  f^^^  is  known  at  the  k  level,  and  (v^.  is  an  added,  fourth 

order  dissipation  term  in  c  as  discussed  in  Beam  and  Warming.^ 
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Typically  X^®,  Y^®  are  approximated  as  follows: 


X 

Y 


C 

C 


o 


o 


^Vl.k  - 

^Vi,k  -  Vi.k)/2 


While  X^®,  Y^®  are  extracted  from  the  nonlinear  differential  equations  (1)  and 

(2). 


This  allows  the  nonlinearity  to  be  maintained.  It  is  important  to  note  that 
the  cell  areas  at  the  k  and  k+1  levels  are  assumed  to  be  known.  This  can  be 
accomplished  in  several  ways.  Steger  and  Chaussee  construct  polar  grids  about 
two  individual  circles  whose  circumference  is  the  arc  length  of  the  body  to  be 
meshed.  The  grids  on  these  two  circles  have  the  same  spacing  in  the  radial 
direction,  (usually  an  exponential  stretching  to  cluster  points  near  the 
body).  However,  the  two  circles  differ  in  grid  spacing  in  the  circumferential 
direction.  One  circle  has  points  placed  in  equal  spaced  increments  around  the 
circumference.  Cell  areas  are  extracted  from  these  polar  grids  such  that  when 
near  the  body,  cell  areas  are  calculated  using  the  nonuniform  circumferential 
distribution  of  points.  When  far  from  the  body,  areas  from  the  uniform 
circumferentially  distributed  circle  are  used.  A  smooth  function  transitions 
from  one  area  type  to  the  other  for  points  between  these  two  extremes. 

The  solution  algorithm  based  on  these  equation's  exhibits  difficulties 
for  geometries  with  slope  discontinuities  and  regions  of  reverse  (concave) 
curvature.  Discontinuities  can  propagate  into  the  grid  interior  with  undesir¬ 
able  results.  Drawing  from  experience  in  other  hyperbolic  systems  it  is 
possible  that  these  problems  may  be  circumvented  by  carefully  including  other 
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forms  of  dissipation  without  drastically  compromising  the  orthogonal  qualities 
of  the  grid.  If  an  analogy  is  made  between  marching  in  time  and  marching  in  n, 
it  seems  reasonable  that  adding  both  temporal  and  spacial  dissipation  would  be 
beneficial  in  insuring  smoothness.  To  accomplish  this  a  more  general  class  of 
integration  in  n  was  chosen: 


^k+l  — 
An 


(l-a) 


/•9Rs 

'W' 


+  a  (i£) 
3n 


'k+1 


When  a  =  1,  the  original  backward  differenced  scheme  is  obtained.  From  a  Taylor 
expansion  it  can  be  shown  that  the  numerical  error  term  for  this  scheme  is: 
ii(l-2o.)  An0  *  higher  order  terms.  Note  that  for  o  M  the  Integration  Is  a 
trapezoidal  type  and  is  fonnally  second  order  in  n.  For  a  >  i  the  error  tem 
has  a  dissipative  effect  in  the  n  direction.  Rewriting  the  algorithm  in  the 
so  called  "delta"  form,  such  that  increments  in  R.  (R^^^  -  are  solved, 
the  algorithm  becomes 

[I  +  o  B  5^]  (R^^j  -  Rj^)  s 

Adding  a  second  order  implicit  dissipation  term  to  this  algorithm  in  delta 
form  serves  to  augment  the  amount  of  explicit  dissipation  that  may  be  added, 
while  not  rormally  degrading  the  accuracy  of  the  method.  Therefore  the  final 
form  of  the  algorithm  is: 


[1 


((Vjdj)  +  o  B-h  Sj]  -  R^)  .  B-1  +  5^(7jaj)2R^ 
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The  marching  algorithm  resulting  from  the  hyperbolic  equations  provides  an 
efficient  means  for  constructing  a  computational  grid.  Figure  2  is  a  flow 
chart  of  the  computer  program  based  on  this  algorithm.  The  program  begins  by 
initializing  the  variables  and  reading  in  the  user  defined  body  coordinates. 

The  marching  procedure  loop  consists  of  four  steps.  First,  the  areas  at  the 
k+1  level  are  approximated.  Then  a  system  of  equations  is  set  up  to  determine 
the  increments  in  X  and  Y,  This  results  in  a  2x2  block  tridiagonal  matrix 
which  is  solved.  The  increments  in  X  and  Y  are  now  known  at  the  k  level  for 
each  j  station,  and  the  grid  is  computed  at  the  k+1  level  using  these  increments. 
The  process  is  repeated  until  the  KMAX  level  is  reached. 

Table  I  is  a  listing  of  the  subroutines  used  in  a  computer  program  that 
utilizes  the  algorithm  described  above.  A  brief  description  of  the  purpose  of 
each  subroutine  is  also  provided. 
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FLOW  DIAGRA/'I 


Body  Defined 


Compute  Metrics 

Compute  Areas  at  k+1 

Set  up  Right  Hand 
Side  b  vector 

Set  up  2X2  Block  Tri- 
diagonal  matrix  M 


2X2  Block  Tridiagonal 
solution 


yJ*K  i  yjjk  +  avJ 

j,k+l  ^j,k  + 

Write  out  current 
k  level 


FIGURE  2 


TABLE  I 


SUBROUTINE  DEFINITIONS 


SUBROUTINE 

CALLED  FROM 

CALLS  TO 

PURPOSE 

MAIN 

- 

INITIA 

STEP 

OUTPUT 

PROGRAM  CONTROL 

INITIA 

MAIN 

BODIS 

READS  INPUT,  DEFINES  INITIAL 

SARC 

EPSIL 

GRID  SPACING 

BODIS 

INITIA 

- 

READS  BODY  COORDINATES 

SARC 

INITIA 

- 

COMPUTES  BODY  PERIMETER 

EPSIL 

INITIA 

- 

DEFINES  DISTRIBUTION 

STEP 

MAIN 

METRIC 

MARCHES  GRID  GENERATOR 

RHS 

FROM  BODY  TO  OUTER 

FILTRY 

BTRI 

BOUNDARY 

METRIC 

STEP 

- 

RHS 

STEP 

- 

CALCULATES  THE  FORCING  FUNCTION 
(RIGHT  HAND  SIDE) 

FILTRY 

STEP 

GMATRX 

FILLS  THE  TRIDIAGONAL  MATRIX 

GMATRX 

FILTRY 

- 

DEFINES  G  MATRIX 

BTRI 

STEP 

LUDEC 

BLOCK  TRIDIAGONAL  SOLUTION 

LUDEC 

BTRIP 

- 

COMPUTES  L-U 

OUTPUT 

MAIN 

- 

WRITES  DESIRED  INFORMATION  TO 
OUTPUT 

BC 

STEP 

- 

ROUTINE  ADDED  TO  C-GRID  TO 

DEFINE  AFT  BOUNDARY  CONDITIONS. 
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III.  PROGRAM  RESULTS 


Two  computer  programs  were  developed  using  the  algorithm  described  in  the 
proceeding  section.  One  program  produces  an  0-type  computational  grid  and  the 
other  provides  C-type  grids.  Results  from  the  0-type  grid  will  be  shown  first, 
followed  by  a  discussion  of  the  modification  necessary  to  produce  a  C-type  grid 
and  some  example  results.  A  step  by  step  examination  of  the  effects  of  the 
various  input  parameter  options  is  provided  in  the  next  section. 

Figure  3  shows  a  far-field  and  near-field  view  of  one  of  the  first 
attempts  at  fitting  an  0-grid  around  a  NLR  7301  airfoil  section.^  This  shape 
represents  a  relatively  simple  test  condition  except  for  the  region  of  the 
trailing  edge.  Careful  selection  of  the  points  used  to  define  the  geometry  at 
the  trailing  edge  can  significantly  improve  the  results.  Figure  4  is  an 
example  of  an  0-type  grid  around  a  square.  Figure  4a  has  equally  distributed 
body  geometry  points,  whereas,  Figure  4b  shows  the  results  of  clustered  and 
stretched  body  definition  points.  Figures  5  and  6  indicate  the  0-type  grid 
results  for  a  complicated  geometry.  The  geometries  shown  provide  a  very 
severe  test  for  any  grid  generation  scheme  and  is  particulary  difficult  for  a 
marching  (hyperbolic)  procedure. 

0-type  computational  grids  are  widely  used  for  a  variety  of  solution 
procedures.  For  some  applications,  however,  a  C-type  grid  is  preferred.  For 
example,  an  attempt  to  determine  the  viscous  flow  characteristics  around  an 
airfoil  section  and  the  resulting  wake  area  behind  the  airfoil,  require  fine 
grid  spacings  close  to  the  body  surface  and  in  the  wake.  A  C-type  grid  can 
easily  provide  this  kind  of  grid  distribution. 

The  most  obvious  change  required  to  convert  from  an  0-grid  to  a  C-grid  is 
the  establishment  of  a  downstream  boundary  condition.  The  boundary  conditions 
for  a  C-grid,  however,  are  only  slightly  more  difficult  than  for  an  0-grid.  A 
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0-GRlD  ON  NLR730I  AIRFOIL 


FIGURE  3 


12 


0-GRIO^ON  A  SQUARE 


FIGURE  4 
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FIGURE  5 
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separate  subroutine  (BC)  was  provided  to  compute  the  downstream  boundary 
conditions,  but  essentially  X  is  fixed  at  XMAX  and  Y  is  incremented  as  ±  aR. 

The  C-type  grid  algorithm  alters  the  form  of  the  tridiagonal  matrix 
slightly  from  the  form  required  by  the  0-grid  procedure.  Figure  7  shows  the 
general  arrangement  of  the  tridiagonal  matrix  used.  Figure  8  shows  a  typical 
C-grid  for  the  NLR  7301  airfoil.  Many  examples  of  a  C-type  grid  around  the 
NLR  7301  airfoil  are  provided  in  the  next  section  which  includes  a  discussion 
of  the  effects  of  the  various  user  definable  parameters. 
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TRIDIAGIONAL  EQUATION 


TYPICAL  C-6RID  ON  NLR  7301  AIRFOIL 


FIGURE  8 
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IV.  EFFECT  OF  USER  DEFINABLE  PARAMETERS 
It  is  appropriate  to  begin  this  section  with  a  discussion  of  how  the 
"body"  is  defined  for  the  grid  generating  routines.  In  addition  to  the  actual 
body,  as  defined  for  the  0-grid  routine,  a  zero  thickness  extension  of  the 
body  from  the  trailing  edge  to  the  downstream  boundary  is  required  for  the 
C-grid  routine.  The  complete  "body"  (actual  body  plus  zero  thickness  extension) 
is  input  in  the  same  way  as  the  real  body  was  input  to  the  0-grid  procedure. 

That  is,  the  body  geometry  is  input  as  JMAX  Cartesian  coordinates,  starting  at 
XMAX  and  proceeding  around  the  body  in  a  clockwise  direction.  A  more  aesthet¬ 
ically  pleasing  grid  is  obtained  if  the  upper  and  lower  (X,  Y)  locations  for 
the  wake  extension  are  equal. 

This  section  includes  a  description  and  discussion  of  the  effects  of 
varying  the  user  defined  input  variables.  It  is  necessary  to  define  a  baseline 
condition  from  which  variations  of  individual  parameters  will  be  discussed. 

Table  II  defines  the  baseline  parameters  used  to  determine  a  grid  for  the 
NLR  7301  airfoil  section  which  will  be  used  for  most  of  the  examples  of  this 
section.  Figure  9  shows  a  far  field  and  a  near  field  view  of  the  grid  produced 
using  these  baseline  parameters.  The  variables  in  Table  II  will  each  be 
discussed  in  turn. 

JMAX  defines  the  total  number  of  points  used  to  describe  the  body. 

Probably  more  important  than  the  total  number  of  points  is  the  distribution  of 
these  points.  In  a  practical  situation  both  the  maximum  number  (JMAX)  and  the 
distribution  will  be  governed  by  the  body  geometry  and  the  finite  difference 
program  that  will  be  using  the  grid.  Figure  10  shows  the  baseline  configuration 
(JMAX  =  100),  and  grids  with  25%  fewer  (JMAX  =  75)  and  50%  more  (JMAX  =  150) 
points.  The  distribution  of  the  points  around  the  body  is  approximately 
proportional  in  all  three  cases. 
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TABLE  II 


BASELINE  PARAMETER  SETTINGS 


PARAMETER 

VALUE 

REASON 

JMAX 

100 

Number  of  coordinates  for  body. 

KMAX 

40 

Typical  value. 

DSETA 

.004 

Provides  approximately  5  grid  lines  in  a 
turbulent  boundary  layer  with  RL  =  10 
million. 

SETMX 

6 

Typical  value. 

ESCAL 

.005 

Nominal  value. 

SMU 

.1 

Nominal  value. 

SMDIM 

.5 

Nominal  value. 

ALPHA 

1.0 

Implicit  finite  difference  of  order  1. 

20 


BASELINE  GRID 


FIGURE  9 
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KMAX  defines  the  maximum  number  of  grid  levels  or  contours.  Note  that 
the  body  geometry  definition  constitutes  the  first  physical  contour  (K  =  1). 
Each  succeeding  contour  is  determined  by  marching  radially  from  the  previous 
contour.  The  initial  cell  height  in  the  radial  (n)  direction  is  also  approxi¬ 
mately  specified  by  the  user  (DSETA).  The  program  uses  a  Newton-Raphson 
routine  to  determine  the  values  of  an  exponential  stretching  function  (c)  that 
will  smoothly  vary  from  the  initial  increment  (DSETA)  at  the  body  to  the 
location  of  the  outer  boundary  (SETAMX)  in  KMAX  steps: 

R,^  =  +  DSETA  (1+ 

Figure  11  shows  a  typical  variation  of  A  R  with  distance  from  the  body. 
With  JMAX,  SETAMX  and  DSETA  fixed,  the  aspect  ratio  (length  to  width)  of  the 
grid  cells  away  from  the  body  are  governed  by  KMAX.  Figure  12  shows  the 
effect  of  reducing  the  value  of  KMAX  by  50%  (KMAX  =  20)  and  of  increasing  the 
value  of  KMAX  by  50%  (KMAX  =  60)  over  the  baseline  value.  The  aspect  ratio  of 
the  cells  and  the  rate  of  transition  to  equal  area  cells  are  both  significantly 
influenced  by  the  value  of  KMAX. 

DSETA  defines  the  initial  cell  height  normal  to  the  body  surface.  The 
proper  magnitude  for  DSETA  depends  on  the  purpose  for  which  the  grid  is  being 
created.  For  example,  if  a  detailed  analysis  of  the  flow  conditions  in  the 
boundary  layer  of  a  high  Reynolds  number  flow  is  desired,  several  grid  levels 
will  be  needed  within  this  thin  layer  near  the  body.  On  the  other  hand,  if 
the  flow  is  assumed  to  be  inviscid  (no  boundary  layer)  the  grid  distribution 
near  the  body  may  be  governed  by  other  factors,  such  as  cell  aspect  ratio.  If 
the  total  number  of  contours  (KMAX)  and  the  outer  boundary  location  (SETAMX) 
are  fixed,  then  obviously  having  more  contours  near  the  body  requires  a  wider 
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VARIATION  OF  KMAX 


FIGURE  12 
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KMAX  =  60 


spacing  away  from  the  body.  Figure  13  shows  the  variation  in  grid  spacing  for 
DSETA  =  .001  and  .01  compared  to  the  baseline  value  of  .004. 

The  next  parameter  in  the  list,  SETAMX,  defines  the  approximate  location 
of  the  outer  boundary  for  the  grid.  For  the  baseline  airfoil  it  defines 
(approximately)  the  number  of  chord  lengths  away  from  the  body  for  the  outer 
most  grid  level  (KMAX).  Figure  14  shows  a  comparison  of  SETAMX  =  3,  6  and  10 
with  all  other  parameters  at  baseline  values.  Note  that  the  downstream 
boundary  has  been  maintained  at  six  chord  lengths  for  this  comparison.  A 
practical  grid  would  probably  have  the  downstream  boundary  approximately  the 
same  distance  from  the  body  as  the  other  far  field  boundaries. 

ESCAL  defines  the  rate  at  which  the  grid  distribution  around  the  body 
will  transition  from  the  input  coordinate  distribution  to  a  grid  with  equal 
cell  areas  for  a  given  radial  location.  The  scale  factor  used  to  transition 
to  equal  cell  areas  is  of  the  form: 

Scalej^  *  (1.  -  ESCAL)*^"^  k  =  2,  3,  4 . KMAX 

Figure  15  shows  how  the  value  of  ESCAL  influences  this  scale  factor.  ESCAL  = 

0  means  there  v/ill  be  no  tendency  for  the  grid  cells  to  seek  equal  areas  with 
increasing  distance  from  the  body.  In  this  case  the  radial  lines  tend  to 
march  out  from  the  body  with  the  same  distribution  as  the  input  data;  however, 
the  procedure  will  not  allow  radial  lines  to  coalesce  or  cross.  Figure  16 
compares  results  for  ESCAL  =  0.  and  0.025.  Notice  that  the  ESCAL  =  0.025  grid 
has  very  undesirable  characteristics  along  a  line  of  "disturbance"  from  the 
trailing  edge.  Additional  problems  results  when  the  disturbance  reaches  the 
rear  boundary  of  the  C-grid  routine.  The  bad  grid  points  occur  between  regions 
of  relatively  dense  grid  spacing  and  regions  of  sparse  grid  spacing.  Some 
relief  for  this  problem  can  be  obtained  with  the  following  smoothing  parameters. 
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VARIATION  OF  DSETA 


FIGURE  13 
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VARIATION  OF  SETAMX 


FIGURE  14 


SETAMX  =  10 


EFFECT  OF  ESC 


Ul  U  <  _|  Ld 


FIGURE  15 


VARIATION  OF  ESCAL 


ESCAL  =0.0 


ESCAL  =0.025 


FIGURE  16 
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The  next  input  parameter  for  the  routine  is  an  explicit  smoothing  parameter 
(SMU).  This  parameter  defines  the  amount  of  fourth  order  dissipation  (smoothing) 
to  be  used  to  damp  numerical  oscillations.  An  explanation  of  why  the 

smoothing  is  needed  and  how  it  is  implemented  is  given  in  Section  II  and 
Reference  (1). 

SMU  has  an  adverse  effect  on  orthogonality  and  too  large  a  value  for  SMU 
could  cause  the  numerical  procedure  to  become  unstable.  The  value  of  SMU 
input  represents  the  maximum  amount  of  smoothing  for  any  grid  level.  The 
actual  value  used  in  the  program  varies  from  zero  at  the  body  to  this  maximum 

value  far  from  the  body  at  a  rate  equal  to  the  rate  of  transition  to  equal 
cell  areas. 

=  ^^^INPUT  -  SCALE^) 

For  the  baseline  parameters,  variations  of  SMU  seem  to  have  very  little 
effect  on  the  grid  produced  until  values  of  approximately  4.5  are  input. 

Figure  17  demonstrates  this  by  showing  grids  for  SMU  =  0.0,  0.5  and  4.5.  At 
SMU  =  4.5,  very  obvious  instabilities  are  occurring  at  grid  points  far  from 
the  body.  It  is  important  to  remember  that  the  value  for  ESCAL,  SMUIM  and  the 

input  geometry  definition  can  all  have  a  strong  influence  on  how  SMU  affects 
the  grid. 

SMUIM  is  another  smoothing  parameter.  This  implicit  smoothing  is  not  as 
effective  as  the  explicit  smoothing  (SMU)  but  it  will  not  cause  the  numerical 
instability  that  large  values  of  SMU  does.  Implicit  smoothing  has  the  dual 
effect  of  adding  some  higher  order  smoothing  itself  plus  increasing  the  amount 
of  explicit  smoothing  that  can  be  added  before  the  procedure  becomes  unstable. 
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VARIATION  OF  SMU 


FIGURE  17 


Although  there  is  no  theoretical  stability  limit  to  the  amount  of  implicit 
smoothing  (SMUIM)  that  can  be  added,  there  is  a  practical  limit  because  both 
orthogonality  and  accuracy  of  the  grid  are  affected.  Figure  18  shows  the 
effect  of  SMUIM  on  orthogonality.  It  is  obvious  from  Figure  18  that  large 

values  of  SMUIM  have  a  significant  effect  on  the  propagation  of  disturbances 
throughout  the  grid. 

The  ALPHA  parameter  controls  the  nature  of  the  finite  difference  marching 
algorithm  used  to  march  the  grid  from  the  body  to  the  outer  boundary  (see 
Section  II).  A  value  of  ALPHA  greater  than  1  tends  to  weight  the  procedure  in 
favor  of  the  implicit  method  and  has  the  result  of  improving  the  smoothness  of 
the  grid.  The  trade-off  is  between  the  accuracy  (and  improved  orthogonality) 
provided  by  the  ALPHA  =  i  choice  and  the  "robustness"  of  the  ALPHA  greater  than 
1  choice.  The  grid  characteristics  required  by  the  finite  difference  program 
for  which  the  grid  is  intended  and  the  shape  of  the  body  about  which  the  grid 
is  produced  will  probably  determine  what  value  of  ALPHA  is  appropriate.  The 
payoff  from  using  large,  values  of  ALPHA  are  most  evident  on  complex  bodies. 
Figure  19  shows  the  effect  of  ALPHA  on  the  grid  generated  about  sections  of 
the  X-24C.  For  the  baseline  parameters  on  the  airfoil  section,  the  effect  of 
changing  ALPHA  from  0.5  to  2.0  is  quite  small  (see  Figure  20).  There  is  some 
stretching  of  the  grid  in  the  radial  direction  as  ALPHA  gets  larger.  Very 
large  values  of  ALPHA  (up  to  10)  produce  significant  stretching  in  the  radial 
direction  and  causes  some  loss  of  orthogonality  near  the  body. 
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VARIATION  OF  SMUIM 


FIGURE  18 
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FIGURE  I-  VARIATION  OF  ALPHA  ON  X-24C 
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VARIATION  OF  ALPHA 


ALPHA  sO.5 


ALPHA  =2 


FIGURE  20 
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V.  CONCLUSIONS 

In  conclusion,  the  hyperbolic  grid  generating  procedure  described  is  very 
fast  and  provides  a  good  grid  pattern  with  a  lot  of  user  control.  The  routine 
is  tolerant  to  variations  of  all  input  parameters  ("robust")  except  for  the 
parameter  that  controls  transition  to  equal  cell  areas  (ESCAL)  when  rapid 
changes  in  grid  spacings  around  the  body  are  used.  Tolerance  to  ESCAL  could 
be  improved  by  changing  some  of  the  other  parameters  (SMUIM  for  example) 
and/or  adding  or  deleting  input  points.  Additional  work  is  needed  in  reducing 
the  effect  of  the  disturbance  between  high  and  low  density  grid  regions. 
Another  desirable  feature  for  clustered  grids  would  be  to  delay  the  transition 
to  equal  cell  areas  in  the  clustered  regions  until  about  a  body  length  away 
from  the  body  (n  =  1). 
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APPENDIX  A  -  SMOOTHING  INPUT  COORDINATES 

Most  finite  difference  procedures  for  which  a  grid  is  required,  need  the 
distribution  of  points  around  the  body  to  be  second  order  smooth. 2  In  other 
words,  a  plot  of  grid  point  number  (J)  versus  body  perimeter  length  (S)  should 
be  such  that  the  rate  of  change  of  curvature  is  mathematically  smooth. 
Smoothing  of  the  geometry  coordinates  (which  constitutes  the  grid  distribution 
at  the  surface)  must  be  accomplished  by  a  separate  program  before  the  geometry 
is  input  to  the  grid  generating  procedure.  Any  clustering  of  grid  points, 
such  as  around  the  leading  edge  or  at  shock  locations,  must  also  be  included 
in  the  body  coordinates  read  by  the  grid  generating  routine.  These  clustered 
coordinates  should  also  be  second  order  smooth. 

A  program  called  SMOOTH  was  written  that  allows  the  distribution  of 
coordinates  around  the  body  and  wake  to  be  broken  into  any  number  of  intervals. 
Within  an  interval  any  number  of  points  can  be  defined  such  that  the  distributi 
within  the  interval  and  across  the  end  points  is  smooth.  An  exponential  point 
distribution  is  used  within  the  intervals  and  a  tension  spline  is  used  to 
interpolate  between  the  original  body  coordinates  to  define  the  new,  smooth 
and  (if  desired)  clustered  points. 

The  end  points  of  the  intervals  are  identified  by  grid  point  numbers  (J). 
The  first  interval  always  begins  at  the  lower  surface  at  X  =  XMAX  (J  =  1). 

The  last  interval  ends  at  the  upper  surface  at  X  =  XMAX  (J  =  XMAX).  The 
program  must  also  have  an  interval  defined  at  the  leading  edge,  X  =  XMIN, 
location  to  avoid  the  possibility  of  having  more  than  one  Y  for  a  given  X  in 
an  interval.  The  curve  fitting  routines  require  monotonically  increasing  or 
decreasing  values  for  X  in  a  given  interval.  Other  points  where  intervals 
should  be  defined  are  the  lower  and  upper  locations  of  the  trailing  edge  of 
the  body.  This  would  minimize  the  possibility  of  incorrectly  defining  the 
geometry  in  that  region. 


In  addition  to  defining  the  location  of  the  start  (and  thus  stop)  points 
of  each  interval,  the  user  must  define  the  number  of  points  in  each  interval, 
and  the  grid  spacing  at  each  interval  point.  A  typical  set  of  input  values 
for  the  SMOOTH  routine  is  shown  in  Table  III.  This  Table  is  for  the  NLR  7301 
airfoil  and  wake  with  100  points  in  both  the  original  input  geometry  and  the 
new  smooth  coordinates.  Figure  21  shows  a  plot  of  coordinate  point  (J)  versus 
body  length  (S)  for  the  original  coordinates  and  the  new  smooth  coordinates. 
Table  IV  is  a  listing  of  the  original  coordinates  and  the  smooth  coordinates. 

As  mentioned  above,  a  tension  spline  routine  is  used  to  determine  the  new 
coordinate  values.  The  value  of  the  tension  parameter,  SIGMA,  must  be  provided 
by  the  user.  The  range  of  values  is  from  zero  (normal  cubic  spline)  to 
infinity  (straight  line  between  points).  For  Figure  21  a  value  of  1000  was 
used.  Too  small  a  value  for  SIGMA  can  produce  undesirable  (not  smooth) 
perturbations  in  the  grid  distribution;  whereas,* too  large  a  value  may  produce 
"straight  line"  geometry  definition.  Figure  22  shows  the  coordinate  distribution 
obtained  by  using  the  parameters  in  Table  III  and  SIGMA  =  1000.  Figure  23 
shows  a  coordinate  distribution  and  resulting  grid  for  clustered  points 
between  50  and  60%  chord  on  the  upper  surface.  This  is  typical  of  what  might 
be  desired  to  provide  a  close  grid  spacing  in  the  region  of  a  shock  wave.  The 
value  of  the  parameter  ESCAL  in  the  grid  procedure  will  need  to  be  kept  small 
to  keep  the  grid  points  clustered  at  points  away  from  the  body. 
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TABLE  III 


INPUT  PARAMETERS  FOR  SMOOTH 


STARTING 

J 

NUMBER  OF  POINTS 
IN  INTERVAL 

GRID  SPACING  AT 
POINT  J 

COMMENTS 

1 

16 

0.5 

Lower  surface  rear  boundary  to 
lower  surface  trailing  edge. 

9 

35 

.04 

Lower  surface  T.E.  to  L.E. 

48 

36 

.006 

L.E.  to  upper  surface  T.E, 

92 

16 

.04 

Upper  surface  T.E.  to  rear 
boundary. 

100 

0 

0.5 

Number  of  points  is  not  needed 
here. 

Note:  The  number  of  new  coordinates  defined  is  the  sum  of  the  number  of  points 
in  each  internal  minus  the  number  of  internals,  plus  two;  i.e.,  (16  +  35 
+  36  +  16)  -  5  +  2  =  100. 
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SMOOTH  COOROINATE  DISTRIBUTION 
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FIGURE  21 
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TABLE  IV 


ORIGINAL 


6.flonnn 

o.noooq 

5.46707 

-0.00000 

_■  4.774.^4 

3.87389 

-0.00000 

2.64807 

0.00000 

1.49351  - 


1.23904  0.00000 


1 .09032 


l.OOOCO  -0.00040 


0.91671  0.00186 

-Q-.afiL5a5 _ =J3..  0C17C  I 

0.81745  -0.00909  I 

■0.761P4  .n.  01979 _ ! 

0.69984  -C. 03296 ! 


0.57206  -0.05974  i 

■Q-i»-SlS4S _ =J)_.  C.6_34i_] 

0.46359  -0.07352  j 

J_._415e6  -0 .07610  ' 

0.37230  -0.07745  : 

■132  47  -C.0  77P  1 _ | 

0.29604  -C. 07727  | 

0.26280  -0.07594  ! 


0.23246  -0.07412  : 

0.«_2.0.478 _ .- 0_t0  719 1_' 

0.17945  -0.06955 

0 .156 43  -0.06709 

0.13546  -0.06*61 

0.11634 _ -0.06184 


0.06830 

-0.05257  ■ 

0.05546 

•0.049C5 

0.04369 

-0.0*527  : 

0.03321 

-0.04094  * 

0 .02388 

-0.03620 

0.01606 

_-0. 03090  ■ 

6  .00969 

-0.02553 

0.00497 

-0.01899  ■ 

0.00120 

0.01310 

O.OOISO 

0.01580 

0  .00290 

0.0198C 

5  0  .00570 

0.02750 

^  C. 00830 

0.0326C 

0 .00539 


01 
0.1 

0.02676  i 
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TABLE  IV  CONCLUDED 


ORIGINAL 


SMOOTH 


0.01660 


0.02420 


0.03120 


.04740 


.08100 


0  .12270 


0  .18320 


0.22960 


0.59950 


IVJ 


20 

50 


0  .80400 
C. 84070 


0  .88060 
0  .92400 


0  .96370 
0  .98310 


1.12500 
1  .50000 


5.00000 


0.04380 


0.05020 


0.05440 


.06090 


.06910 


0.07560 


8200 


0.08560 


887 


0.03134 


4169 


1.05447 


0.05493 

0.06306 

n. 0^844 

0.08365 

0.06958 

aS 


.11868 


flilricLI 


07504 


e, 


33933 


0.09034 


42309  0.08987 


0.52223  0.08718 

n  .57843 _ D...082X2 

0.63974  0.07498 


0.76649  0.05342 


0.96047 


Itllllltlll 


0.04540 

0.0370C 


0.00850 

0.00420 


0 


1.04079 

1.09034 


1 

3 


1.34922 

1.49350 


0.00000 
. QQOOO 
0 
0 


2.23817 

2.64805 


>0.00001 

•O.OGOOO 


0 

0 


•0.0  0003 
POO 


•0.0  0000 
•0.0  000.0 


4.77432 

5.46705 

'6.00000 


•0.0  0000 

0.00000 
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